Coarsening Dynamics of Granular Heaplets in Tapped Granular Layers 
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A semi-continuum model is introduced to study the dynamics of the formation of granular heaplets 
in tapped granular layers. By taking into account the energy dissipation of collisions and screening 
effects due to avalanches, this model is able to reproduce qualitatively the pattern of these heaplets. 
Our simulations show that the granular heaplets are characterised by an effective surface tension 
which depends on the magnitude of the tapping intensity. Also, we observe that there is a coarsening 
effect in that the average size of the heaplets, V grows as the number of taps k increases. The growth 
law at intermediate times can be fitted by a scaling function V ~ k z but the range of validity of 
the power law is limited by size effects. The growth exponent z appears to diverge as the tapping 
intensity is increased. 

PACS numbers: 61.43.Gt,45.70.Qj,83.80.Fg 



In recent years, there has been increasing interest in 
the collective dynamics of granular materials. The dis- 
sipative nature of granular materials gives rise to prop- 
erties distinct from those of solids and liquids. Many 
experimental and theoretical attempts have been made 
to seek a fundamental understanding of this granular 
state and especially pattern formation in driven granular 
layers. The experiments include vertically vibrated sys- 
tems |y, , tapped or blown thin films of powders j|, [| , 
and electrostatically driven granular layers || . Al- 
though many papers have been published to explain the 
experiments on tapped granular layers, most of them deal 
with either compactification of thick granular layers @, |) , 
or static properties of granular heaplets on a thin gran- 
ular layer || . The dynamical aspect of the heaping pro- 
cess in a tapped thin granular layer is still not well un- 
derstood. This paper is concerned with the coarsening 
dynamics of granular heaplets in tapped granular layers 
from a theoretical point of view. 

This paper is organised as follows. First we introduce a 
simple model to study this fascinating system by consid- 
ering the energy dissipation and screening effects during 
the tapping process. Despite the simplicity, this model 
is capable of capturing the essential phenomenology, re- 
producing various morphologies of the coarsening pat- 
tern, and showing the way in which heaplets merge as 
observed in experiments |3| . The model is then studied 
numerically and the results indicate some of the relation- 
ships between tapping intensity and the effective surface 
tension of granular layers. Finally analysis of the results 
shows that the average size of the the heaplets V grows 
as a power law with the number of taps k, V ~ k z for 
limited range of k values. The exponent z appears to 
diverge as the tapping intensity is increased. 

Model - Consider a thin layer of N ^S> 1 granular parti- 
cles spread out over a fiat plate. The plate is then tapped 
repeatedly at a low pace with constant shock intensity. 
In one complete tapping cycle, there are two different 
processes each of which requires a different description. 
In the tapping phase, the granular layer is perturbed by 



an external shock and granular particles acquire kinetic 
energy to allow them explore the phase space and hop 
on to neighbouring sites. The hopping range of the gran- 
ular particles depends on the amount of kinetic energy 
received by each of the individual particles. The amount 
of energy supplied to the system is controlled by the di- 
mensionless acceleration V = U/g At, where U is the ve- 
locity of the plate when it is in motion and At is the 
time interval over which it is in motion. Note that U 
is also the vertical velocity gained by an elastic parti- 
cle sitting on the plate when it is tapped. Our particles 
are not however elastic. Therefore, the vertical velocity 
gain is all, where < a < 1 is the departure coefficient 
analogous to the coefficient of restitution and it charac- 
terises the degree of dissipation of the system. After the 
tap, the particles undergo a ballistic flight and fall back 
again onto the static plate where they relax subject to 
avalanche dynamics before the next tap. In this relax- 
ation phase, the granular particles may stay immobile or 
move about depending on the local slope. If the local 
slope is less than a critical slope then the profile remains 
stationary. Conversely, if the local slope is greater than 
the critical slope, matter moves down the slope collec- 
tively as an avalanche, until the slope is less than the 
critical slope. 

Let the area density n(x, k) be the number of granular 
particles above a unit area of the plate after k taps. Here 
x = (x,y), where x and y are orthogonal coordinates 
parallel to the plate. The average diameter, D of the 
granular particles is chosen to be 1, so that n(x, k) now is 
equal to the local thickness of the granular layer, so long 
as there is no compactification throughout the tapping 
process. Of course this is just an idealisation, the packing 
fraction of the granular material will be changed |J if the 
granular layer is thick or it is intentionally prepared in a 
low-density state, but such cases will not be considered 
here. 

We can approximate the acceleration A p of the plate's 
movement due to tapping as a sum of M delta functions 
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Ap(t) = Y,U8(t- (Jfe-1) 



(1) 



fe=i 



where the velocity U gives the intensity of the shocks 
which are assumed uniform over the whole plate, r spec- 
ifies the time interval between taps. Here r must be 
greater than any time scale in the relaxation process, this 
is to ensure that next tap only occurs after the system is 
fully relaxed. In the dilute limit, a single inelastic grain 
sitting on the plate experiences a net force 
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where a is the departure coefficient mentioned earlier. 
During a small interval [tk — At/2, tk + Af/2], where tk 
is the short hand for tk = (k — l)r, the time of fc-th tap, 
and the maximum velocity transferred to the particle is 



Av z = 



t k +- 



dt = all - g At. 



(3) 



In order for the particle to hop, Av z must be positive, or 
all > Vq = g At so that external acceleration overcomes 
the gravitational force. In what follows this is always the 
case. 

In the continuum limit one has to make two crucial 
modifications, a must be replaced by an effective de- 
parture coefficient of the granular bulk, a is expected 
to be a monotonic decreasing function of n(x, k). This 
is because when the number density is high, inter-grain 
collisions occur more often and more energy is dissipated, 
so that grains depart with a smaller departure velocity. 
While the precise form of a(n) is unknown, in this paper 
it is taken to be 



a(n) 



H^TT, n > 1 

ocqu^ n < 1, 



(4) 



for reasons of simplicity and to make comparison with 
Duran's model ||. Here, n is the average density of the 
system and ao is the effective departure coefficient for 
average density n. Also, the second term in Eq.(|^) needs 
to be changed. In the continuum limit, according to Du- 
ran Q there is an effect dependent on the position of a 
grain in the heap. If the grain is not at the top of the heap 
it supports a fraction of the weight of the grains above 
it, hence its effective mass is increased. As a result, in 
order for a granular particle sitting on the inclined side 
of the granular layer to hop, it requires a larger velocity 
kick than a particle sitting on the flat region of the layer. 
This screening effect can be incorporated into Eq.(j^) by 
replacing m in the second term with an effective mass 



m* = {{iit — n)p sin C /D) m. 



(5) 



Here D is the average grain diameter and it is set to I 
henceforth, n-r is the altitude of the nearest peak in the 



corrugated granular layer, p is a parameter of unknown 
value which is set equal to 5 in reference Q] and this value 
appears to give a match to experimental results. 9 C is the 
angle of repose and is close to the value of it/ 6 M. After 
these modifications, Eq.(||) now can be written as 

Ar. ( A - t a , i, ),, >UI f) a .. | 
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(jit — n)p sin 9 C n > I 
(tit — n)p sin# c n < 1, 



where A 



annlJ 



is the tapping intensity. 



There are close similarities between our model and Du- 
ran's H. In both models, the resulting pattern formation 
is due to the competition between the upward hopping 
motion of particles and the downward screening effect 
of avalanches. However, the mechanism causing upward 
motion is different in the two models. Duran [|3j and 
Shinbrot [^0) conjectured that the hopping of granular 
particles is due to the upcoming air-flux through the 
porous bed of the granular layer. The velocity of the 
air-flux at height n can be approximated by Darcy's law 
v a ir = Kp/n, where K is the permeability of the granular 
bed, p is the pressure different across the granular bed. 
The velocity of the air-flux v a i r must be greater than 
the free fall velocity of a granular particle Vf in order 
to eject this particle. However, in our model we adopt 
a simpler picture where kinetic energy is transferred to 
the granular particles via direct collisions between plate 
and particle and between particle and particle. Energy 
is dissipated via the effective departure coefficient. The 
form of the effective departure coefficient is chosen as 
in Eq.(||) so that the velocity transfer from plate to the 
particle has the same 1/n dependency (n > I) as v a i r in 
Darcy's law. Obviously, we have the freedom to choose 
other functional forms of a(n), the simple choice here 
is just to achieve comparison with Duran's model. As 
observed from simulation, however, it turned out that 
the exact form of a(n) is unimportant so long as it is a 
monotonic decreasing function. 

Now it remains to specify the hopping process in the 
tapping phase. When the granular layer is tapped, at 
an unstable site where Av z > 0, particles can hop to 
any neighbouring site within a circle defined by a max- 
imum horizontal hopping range, R max {Av). Here we 
assume on the unstable site, all the granular particles 
are re-distributed. Each grain hops with equal prob- 
abilities in any direction and with a random hopping 
range, R < R ma x- We can estimate the maximum hop- 
ping range R max (Av) by the following simple arguments. 
The maximum velocity received by a particle during one 
tap is Av z vertical to the plate. In this case R max {Av) 
is zero for the motion is strictly vertical. It is unlikely 
that particles will always hop in the vertical direction, 
an inter-particle collision can change the hopping direc- 
tion to any angle. Assume that particle is ejected with a 
speed Av z cos <f> with <fi denoting the angle from the verti- 
cal axis of the plate. Then, the horizontal hopping range 



R = Av z cos cj> sin <f> t' 
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<Av z )\ (7) 
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where t' is the flight time of the particle in vacuum. R 



is maximum when <f> 

for Rmax 

is given by 



7r/6, so that the estimated value 



Rmax {Av z 
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where 
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The estimated i? m0 a; is very crude and does not take 
into account the details of the air-particle and particle- 
particle interactions except phenomenologically. How- 
ever, one can consider £ in R max as an adjustable pa- 
rameter, which varies from one experiment to another 
experiment depending on the roughness of the material 
used in experiment and humidity of the surrounding en- 
vironment. 

There are several ways to describe the avalanche pro- 
cess during the relaxation phase. Here we use a slope de- 
pendent diffusion equation. We start from the equation 
of continuity dtfi = — V • J with the constitutive current 
J given by 



-?7(/3|Vn| 2 - l)Vn, 
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where l/\fff = tan6* c is the critical slope and r) sets the 
diffusion rate. The current is chosen in such a way that 
when the gradient is greater than the critical slope, mass 
current flows down-hill to smooth out density fluctua- 
tions, but no mass is transferred if the gradient is less 
than the critical slope. 

Simulation - The model is studied on a N x N square 
lattice with a periodic boundary condition. Initially, the 
granular layer is prepared by assigning a random height 
between [0, 5] to each site and letting the granular layer 
relax. At the beginning, the fluctuation in n(x, 0) is con- 
siderably smaller then the height of the heaplets formed 
later. On each tapping cycle, the velocity of particles at 
each site is calculated and the density number is updated 
according to rules defined by Eq.@ and Eq.(|J). After 
which the equation of continuity and the corresponding 
constitutive current equation Eq.(|l0|) are solved and it- 
erated numerically until there is no more mass is trans- 
ferred at each site. 

The following density plots show typical results from 
the simulation. Fig|l| shows the density plots of n(x, k) 
at different values of k, the number of taps, for two set 
of parameters. The top panel shows extended patterns 
of ridges and corresponds to a smaller tapping inten- 
sity (A — 6.0), while the lower panel shows localised 
heaplets and corresponds to a greater tapping intensity 
(A = 16.8). Both simulations coarsen when k increases, 
eventually reaching their final stage where the entire sys- 
tem is a single granular heap. 

One can see from Figjl] that as the tapping intensity is 
increased, the patterns favour isolated circular heaplets 



which suggests that an effective surface tension can be 
defined which increases with the tapping intensity. Since 
there are no attractive forces between granular particles, 
the granular layer cannot have a true surface tension. 
This "surface tension" effect is entirely due to the sys- 
tem trying to maximise the local energy dissipation by 
decreasing the value of the effective departure coefficient 
so that particles are less mobile. In order to maximise 
the local energy losses, particles have to be as close to 
each other as possible, for this will increase the num- 
ber of inelastic collisions between particles. Therefore, 
with the constraint that the local slope should not ex- 
ceeds the critical slope, the global attractor of the pat- 
tern is a single large heap, where local energy loss is the 
largest. Akiyama et al. Clement |fj"T) and Duran || 
have discussed the effects of convection on heaping pro- 
cess of granular layers. Duran Q in particular suggests 
that the "surface tension" is brought about by the con- 
vective forces dragging particles from the surroundings 
into the granular heaplets. However, in our simulation 
it is due to the maximisation of energy losses, since in 
our model the granular layer does not contain any infor- 
mation about the internal interactions of the sand heap, 
therefore convective forces are not taken into account. 




FIG. 1: Two simulations of coarsening dynamics of tapped 
granular layers corresponding to different tapping intensities. 
The simulation images (a)-(d) and (e)-(h) start with the same 
initial configuration. Figures (a)-(d) use a smaller tapping in- 
tensity (A = 6.0, £ = 0.2207), and Figures (e)-(h) correspond 
to a stronger tapping intensity (A = 16.8,^ = 0.2207). (a) 
and (e) are taken at k = 5, (b) and (f) at k = 10, (c) and (g) 
at k = 20, (d) and (h) at k = 30. 




FIG. 2: Merging of two granular heaplets. A = 16.0, £ = 
0.2207. When two heaplets touch each other, smaller heaplet 
is sucked into the larger one. (a) k — 15 two granular heaplets 
meet, (b) k = 30 smaller heaplet is sucked into the larger one. 
(c) k — 50 small heaplet disappears. 
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Nevertheless, the simulation result does show the 
Laplace- Young pressure effect suggested by reference j| . 
A typical example is shown in Fig.||. (The parameters 
here are A — 16.0 and £ = 0.2207, but similar results are 
found over a wide range of values of these parameters.) In 
the center of the figure there are two connected granular 
heaplets. After several taps the smaller heaplet merges 
with the larger heaplet due to matter moving along the 
connecting neck. 

We are also interested in the dynamics of the size of 
the heaplets. The relevant length scale measure I is the 
average thickness of the site hi weighted by the local 
volume, Vi, and the average volume of the heaplet is V — 
I 3 . Since we know that hi = rii and the local volume of 
the site is just Vi — riicr, where a is the area element of 
the lattice site, it follows that 
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FIG. 3: Log-log plot for characteristic heaplet size V against 
number of taps k for varying tapping intensities A. The solid 
line corresponds to A — 18.0, the two dashed lines to A = 17.6 
the dash-dotted line to A = 16.0, and the dotted line to A = 
14.4. The parameter used in the simulation is £ = 0.2207. 
All curves correspond to a lattice size L\ x L\ = 80 x 80 and 
are averaged over ten different runs except for the marked 
dashed curve which corresponds to L2 x L2 = 120 x 120 and 
is averaged over five runs. Each curve saturated after a long 
time, this can be clearly seen from the A = 18.0 solid line. 
The marked curve corresponds to the larger lattice size shows 
a wider range of the linear region and a larger saturation value 
of V. 

There is some limited evidence for power law behaviour 

V ~ k z in the intermediate region of the log-log plots of 

V versus k in Fig.||. There are three regions in each 
plot. The early region is where the granular layer is ran- 
domly distributed, the number density fluctuation is not 
large enough to trigger coalescence. The late-time re- 
gion is where all the heaplets have joined into one big 



heap, leaving some remaining individual particles hop- 
ping randomly about and occasionally encountering the 
large single heap. The intermediate region is also the 
fast-growing region where heaplets grow and merge into 
each other. The range of the intermediate region grows 
and remains linear as the system size increases as can be 
seen from the two curves with A = 17.6. This shows that 
system size limits the range of the intermediate power 
law region. For very large system sizes we would ex- 
pect on this basis that the intermediate power law region 
would have a greater range and remain linear in the log- 
log plot. One can easily estimate the increase in the 
average volume of the heaplets V by the following sim- 
ple argument. As the saturated value of V corresponds 
to a state where a single heap contains approximately 
all the grains in the system, it should be proportional 
to nl/ 2 , here the average number of grains per site n is 
a conserved quantity. Then, increment in V is given by 
A log 10 V — 2A log 10 L. In Fig.| the dashed line corre- 
sponds to a lattice size L\ x L\ = 80 x 80 and the marked- 
dashed line to L2 x L 2 = 120 x 120, and the vertical shift 
in the graphs is expected to be A log 10 V = 0.352 which 
may be compared to 0.360 from the simulation. 
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FIG. 4: Exponent z for varying tapping intensity A. For small 
A (A < 14), z is roughly constant (z ~ 1.0), but appears to 
diverge for A > 14. Each value of z is averaged over ten 
different runs. 

The exponent z in the power law depends on the tap- 
ping intensity A as is shown in Fig.[| As we can see z 
is roughly constant, z ~ 1.0, for low tapping intensity. 
When the tapping intensity is increased A > 14, z sud- 
denly increases, and appears to diverge near A = 18. For 
these large values of the tapping intensity, we observed 
that the system appears to be in a flat configuration for 
some time and suddenly a single heap formed within a 
few taps. For still greater tapping intensities, the heaplet 
patterns disappear and the layer remains flat with a small 
random variation superimposed. This pattern-disorder 
transition occurs at a not very precisely defined critical 
tapping intensity of A ~ 18. (This critical intensity ap- 
pears to depend slightly on the initial conditions which 
is the reason for our statement that it is not precisely de- 
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fined.) Above this critical tapping intensity, no patterns 
are observed. This is due to the fact that the pertur- 
bation is so strong that no site is stable, and each site 
constantly undergoes re-distribution on each tapping cy- 
cle. Although the power law parametrisation is based on 
limited evidence, it does provide a useful parameter z to 
distinguish the pattern- forming region (A < 18) from the 
non-pattern- forming region (A > 18). 

We have introduced a simple model of the heaping pro- 
cess in a tapped granular layer. The model is capable of 
reproducing the essential morphologies of tapped gran- 



ular systems. Qualitatively the effective surface tension 
of the granular heaps is closely related to the tapping 
intensity, and it is shown that there is no need for con- 
vective forces for this effective surface tension to exist. 
The length scale of the system coarsens according to the 
power law I ~ k z in a limited range. The exponent z ap- 
pears to diverge as tapping intensity is increased and pro- 
vides a useful parameter for distinguishing the pattern- 
forming and the non-pattern-forming regions. 

We thank Philip Cheung for helpful discussions. 
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